Here, we use Mathematica to graph the phase portrait of  a linear system which is an exercise in the textbook. Remember that we can solve all linear systems!

({{x'= 3x +4y}, {y'= 4x-3y}})

Do each part separately before combining them together, as to make sure there are no errors.

graph1=VectorPlot[{3x+4y,4x-3y},{x,-5,5},{y,-5,5}]

The direction field is working. Now, try some solution curves.

NDSolve[{x'[t]==3x[t]+4y[t],y'[t]==4x[t]-3y[t],x[0]==1,y[0]==1},{x[t],y[t]},{t,-10,10}]

The numerical solution command is working. Based on the qualitative analysis, choose appropriate initial conditions as to represent all types of solutions.

sol1=NDSolve[{x'[t]==3x[t]+4y[t],y'[t]==4x[t]-3y[t],x[0]==0,y[0]==1},{x[t],y[t]},{t,-10,10}];
graph2=ParametricPlot[Evaluate[{x[t],y[t]}/. sol1],{t,-3,3}, PlotRange->{{-3,3},{-3,3}}]

The solution with initial condition (0,1) is working. Notice we have used both negative and positive time to show the solution in both directions, not just the increasing time. For linear systems, this should always work!
Now, combine them together.

graph1=VectorPlot[{3x+4y,4x-3y},{x,-5,5},{y,-5,5}] ;sol1=NDSolve[{x'[t]==3x[t]+4y[t],y'[t]==4x[t]-3y[t],x[0]==0,y[0]==1},{x[t],y[t]},{t,-10,10}];
graph2=ParametricPlot[Evaluate[{x[t],y[t]}/. sol1],{t,-3,3}, PlotRange->{{-3,3},{-3,3}}];
Show[{graph1,graph2},PlotRange->{{-5,5},{-5,5}},AxesLabel->{"x","y"}, Axes->True]

Now, add more solutions.  Based on the qualitative analysis, choose appropriate initial conditions as to represent all types of solutions.  For example, make sure to pick initial conditions on the straight-line solutions, if any.

graph1=VectorPlot[{3x+4y,4x-3y},{x,-5,5},{y,-5,5}] ;
sol1=NDSolve[{x'[t]==3x[t]+4y[t],y'[t]==4x[t]-3y[t],x[0]==0,y[0]==1},{x[t],y[t]},{t,-10,10}];
graph2=ParametricPlot[Evaluate[{x[t],y[t]}/. sol1],{t,-3,3}, PlotRange->{{-3,3},{-3,3}}];
sol2=NDSolve[{x'[t]==3x[t]+4y[t],y'[t]==4x[t]-3y[t],x[0]==0,y[0]==-1},{x[t],y[t]},{t,-10,10}];
graph3=ParametricPlot[Evaluate[{x[t],y[t]}/. sol2],{t,-3,3}, PlotRange->{{-3,3},{-3,3}}];
Show[{graph1,graph2,graph3},PlotRange->{{-5,5},{-5,5}},AxesLabel->{"x","y"}, Axes->True]

To do this more efficiently, we can define a function of solutions, gensol, with independent variables initial conditions, x0 and y0. Then plug in values for x0 and y0, to get different solutions. Finally, graph all using a Table command. Based on the qualitative analysis, choose appropriate initial conditions as to represent all types of solutions.  For example, make sure to pick initial conditions on the straight-line solutions, if any.

graph1=VectorPlot[{3x+4y,4x-3y},{x,-5,5},{y,-5,5}] ;
gensol[x0_,y0_]:=NDSolve[{x'[t]==3x[t]+4y[t],y'[t]==4x[t]-3y[t],x[0]==x0,y[0]==y0},{x[t],y[t]},{t,-10,10}];
sol[1]=gensol[0,1];sol[2]=gensol[0,2.5];sol[3]=gensol[0,-1];sol[4]=gensol[0,-2.5];sol[5]=gensol[1,0];sol[6]=gensol[3,0];sol[7]=gensol[-1,0];sol[8]=gensol[-3,0];sol[9]=gensol[2,1];sol[10]=gensol[-2,-1];sol[11]=gensol[-1,2];sol[12]=gensol[1,-2];
graph2=ParametricPlot[Evaluate[Table[{x[t],y[t]}/.sol[i],{i,12}]],{t,-3,3}, PlotRange->{{-3,3},{-3,3}}];
Show[{graph1,graph2},PlotRange->{{-5,5},{-5,5}},AxesLabel->{"x","y"}, Axes->True]